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20.  Abstract  (Continued) 

^-latitude-longitude  grid  containing  2J  X  (J  -  1)  +  2  data  points.  Numerical 
results  are  presented  to  demonstrate  the  method's  accuracy  and  efficiency. 
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Double  Fourier  Series  Solution 
of  Poisson's  Equation  on  a  Sphere 


I.  IM'KODI  CTION 

Advances  in  numerical  simulation  and  predict  :•  <n  in  d.sripliucs  as  diverse  as 
geophysical  fluid  dynamics,  heat  transfer,  and  mad'  j.  and  plasma  f'i.vsics  have 
generated,  over  the  years,  considerable  interest  in  H>e  method  i.f  solution  f«>r 
Poisson-type  equations,  Most  research  in  this  area,  however,  rloes  not  deal  with 
spherical  geometries.  In  fact,  of  a  list  of  150  ariirles  on’  Fast  Elliptic  Solvers" 
compiled  by  Schuman,  ^  only  two^  ^  deal  with  the  spherical  geometry.  We  present 

in  this  paper  a  new  numerical  method  for  the  solution  of  Poisson's  equation  on  a 

2  3 

sphere.  In  contrast  to  previous  methods  of  Swarztraubcr  and  Yee  which  are 

based  on  finite-difference,  this  method  is  based  on  truncated  double  Fourier  series 

4  5 

on  spheres  (for  example,  Ors/.ag  ;  Boer  and  Steinberg'  ).  The  use  of  finite  double 
(Received  for  publication  28  October  1980) 

1.  Schumann,  V. ,  (ed.  )  (1978)  Fast  Elliptic  Solver’s,  Proceedings  of  the  GAMM- 

workshop  on  fast  solution  methods  for  the  discretized  Poisson  equation, 
Karlsruhe,  1977.  Advance  Publications,  London  WC2,  England,  289  pp. 

2.  Swa rztraube r,  P.  N,  (1974)  The  direct  solution  of  the  discrete  Poisson  equation 

on  the  surface  of  a  sphere,  J.  C'om put.  Phys.  15:40-54. 

3.  Yee,  S.  Y.  K.  (1970)  An  efficient  method  for  a  finite-difference  solution  of  the 

Poisson  equation  on  the  surface  of  a  sphere,  J,  Cumput,  Phys,  22  :21 5-228. 

4.  Orszag,  S.A.  (1974)  Fourier*  series  on  spheres,  Mon.  Wra.  Rev.  102:50-75. 

5.  Boer,  G.J.,  and  Steinberg,  I,.  (1975)  Fourier  series  on  spheres,  Atmosphere 

b3: 180-  191,  - 
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fracnaw  p*j*  bunk-mot  filmed 


Fourier  series  on  a  rectangle  is  not  new,  but  the  formulation  on  a  sphere  reported 
here  is.  In  particular,  Yee's  interpretation  of  Fourier  series  on  spheres  has 
been  implemented  to  narrow  what  Swarztrauber  referred  to  as  ’’the  theoretical  gap 
which  exists  between  the  states  of  the  art  for  discrete  spectral  approximations  on 
a  sphere  and  on  a  rectangle,  "  and  to  obtain  a  Galerkin  approximation  to  the  solu¬ 
tion  of  Poisson's  equation  on  a  sphere. 

Let  u(fl,  X)  be  a  scalar  function  on  a  sphere  where  the  location  of  a  point  is 

specified  by  co-latitude  0  <  6  <  n  and  longitude  0  <  X  <  2 jt.  It  may  then  be  repre- 

8 

sented  by  double  Fourier  series  of  the  form 


u(0,  X) 


00  oo 

£  X  u£  mI(l  "  s)  cos  to  +  s  sin 

m  =  -oo  t  =0 


(la) 


where 


u 


i ,  m 


IT 


0 


2tt 

J  u{0,  X)  e  dX 
0 


[(1  -  s)  cos  f  9  +  s  sin  £  0]  d0 


(lb) 


j  1  if  (  =  0 
(  2  otherwise. 


and  s  =  0  if  m  is  even,  s  =  1  if  m  is  odd.  Since  the  transform  pair  Eqs.  (la)  and 
(lb)  is  amenable  to  Fast  Fourier  Transform  (FFT),  we  have  at  our  disposal  an 
extremely  efficient  means  of  obtaining  a  numerical  solution  for  Poisson's  equation 
on  the  surface  of  a  sphere. 

The  procedure  involves  the  expansion  of  the  dependent  variables  in  truncated 
double  Fourier  series,  the  substitution  of  the  truncated  series  for  the  dependent 
variables  in  the  Poisson  equation,  the  application  of  the  Galerkin  approximation  to 
obtain  in  Fourier  space  a  number  of  linear  algebraic  systems,  the  solution  of  these 
systems,  and  the  inverse  transform  of  the  solution  in  Fourier  space  back  to  physi¬ 
cal  space. 


6.  Yee,  S.  Y.  K.  (1980a)  Fourier  series  on  spheres-a  geometric  perspective, 

AFGL-TH-80-0021,  Air  Force  Geophysics  Laboratory,  Bedford,  MA.,  16  pp. 

7.  Swarztrauber,  P.  N.  (1979)  On  the  spectral  approximation  of  discrete  scalar 

and  vector  functions  on  the  sphere,  SIAM  J.  Numer.  Anal.  16:934-949. 

8.  Yee,  S.  Y.  K.  (1980b)  Studies  on  Fourier  series  on  spheres,  Mon.  Wea.  Rev. 

108:676-678. 
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2.  ORDINARY  DIFFERENTIAL  EQUATIONS  WITH  VARIABLE 
COEFFICIENTS 


Our  goal  is  to  seek,  using  the  FFT  technique,  an  accurate  numerical  solution 
to  the  Poisson  equation  on  the  surface  of  a  unit  sphere  (0  <  6  ^  ff,  0  ^  A  <  2 ff). 


sin  9  9 6  Sm  0 


au 

90 


2 

sin  0 


3X2 


(2a) 


Here  the  forcing  function  f(R,  A)  must  satisfy  the  compatibility  condition  (for  exam- 
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pie,  Berg  and  McGregor  ) 


/ f(0,  A)  sin  0  ds  = 
S 


0 


(2b) 


over  the  surface  of  the  sphere  S. 

The  first  step  in  our  solution  method  is  the  transformation  of  this  partial  dif¬ 
ferential  equation  (PDE)  to  a  number  of  ordinary  differential  equations  (ODE),  a 
standard  approach  in  the  method  of  solution  for  PDE*  A  brief  outline  of  this  proce¬ 
dure  is  included  here  for  completeness.  We  decompose  u(0,  A)  along  a  given  latitude 
by  a  truncated  Fourier  series  of  the  form 


M 

u(e,A)  =  £  um(o) eim*  » 

m  =  -M 


where 


_c 

K 


K 

X  u(0,  Afc) 
k=l 


e 


-imXk 


(3a) 


c  =  1  for  m  »  0  or  ±M 
c  =  2  otherwise 


(3b) 


are  complex  Fourier  coefficients,  M  =  K/2  is  the  cutoff  wavenumber,  A^  =  2?rk/K, 
and  K  is  the  number  of  data  points  along  a  latitude  circle.  Since  a  latitude  circle 
degenerates  into  a  single  point  at  the  poles,  for  u  (0,  A)  to  have  a  single  value  u(P) 
there,  we  must  have  the  pole  conditions 


9.  Berg,  P.W.,  and  McGregor,  J.  L.  (1966)  Elementary  Partial  Differential 
Equations.  Holden-Day,  San  Francisco,  CA„ 
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(4a) 


!u(P)  for  m  =  0 

0  for  m  *  0  ,  (4b) 

where  P  =  0  or  it,  and  |m|  <  M.  A  Galerkin  approximation  to  Eq.  (2a),  represen¬ 
ted  by  a  set  of  ODE  with  variable  coefficients,  is  obtained  if  we  substitute  for  u 
and  f  in  Eq.  (2a)  by  truncated  series  of  the  form  of  Eq.  (3a),  multiply  the  resulting 
equations  by  exp  (-imA),  and  then  integrate  over  X  for  the  interval  0  <  X  <  2^: 


sin  0  d0 


sin  8 


d_ 

dO 


um<0)  - 


•  2  . 
sin  0 


um<0)  =  fra<«> 


(5) 


where  0  <  S  ^  |ml  <  M.  For  convenience  of  discussion,  we  rewrite  Eq,  (5) 
in  the  form 


(0)  + 


cos  6 
sin  0 


u»  (0)  - 
m 


•  2  a 
sin  tf 


u  (6) 

m 


f  <0) 

m 


(6) 


Equation  (6)  together  with  Eq.  (4)  constitutes  a  set  of  one-dimensional  Dirichlet 
boundary-value  problems  (BVP).  Of  this  set,  2M  equations  have  homogeneous 
boundary  conditions  and  one  equation  has  inhomogeneous  boundary  conditions  given 
by  Eq.  (4a).  Thus  the  problem  of  solving  a  two-dimensional  Poisson  equation  be¬ 
comes  that  of  solving  a  set  of  one-dimensional  variable  coefficient  BVP  of  the 
Helmholtz-type.  The  pole  conditions  may  now  be  considered  as  boundary  condi¬ 
tions.  For  m  =  0,  however,  Eq.  (6)  has  no  unique  solution  because  the  pole  values 
u(P)  in  Eq.  (4a)  are  among  the  unknowns  to  be  sought  in  Eqs.  (2a)  and  (2b).  This, 
of  course,  indicates  nothing  more  than  that  solutions  to  PoissonTs  equation  on  a 
sphere  can  only  be  determined  to  within  an  additive  constant.  We  shall  return  to 
this  point  later. 


3.  SPECTRAL  ALGEBRAIC  EQUATIONS 

The  next  key  step  in  our  solution  method  is  based  on  the  fact  that  u  (0)  is  a 
2^-periodic  even  (odd)  function  for  even  (odd)  m,  and  we  may  approximate  um(^) 
truncated  half-range  expansions  of  cosine  (sine): 
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u  (6)  = 
m 


Z  uf  cos  f  0  for*  even  m 


Eu.  sin  ( 9  for  odd  m  , 

v  m  ' 


where 


-i  0  j  £  u0(fl.)  cos  se-\  +  “^^(0  -  o)  +  (-i)(  u(0  tt)1 


T  1C  u  (0  )  cos  f  0.  for*  even  m  ^  0 
J  m  .1  ,1 


—  u  (0.)  sin  f  0.  for  odd  m 
J  m  j  j 

.i  1 


c  =  1  if  i  =  0  or  L,  c  2  otherwise,  L  =  J  is  a  cutoff,  0.  -  jjr  /  J,  and  (J-l)  is 
the  number  of  data  points  between  poles. 

We  consider  first  the  constraints  on  the  use  of  these  expansions  to  solve 
Eq.  (6).  For  m  =  0,  Eq.  (G)  requires  that  u{j(0)/sin  0  be  finite  for  all  0,  including 
0  0  and  n.  The  cosine  expansions  in  Eq.  (7a)  obviously  satisfy  this  requirement. 

For  m  *  0,  Eq.  (6)  requires  at  least  that  um(0)/sin  0  be  finite  for  all  0  (C)rs?ag4'. 
In  the  case  of  odd  mf  this  condition  is  automatically  satisfied  by  the  sine  expan¬ 
sions  in  Eq.  (7a).  In  the  case  of  even  m  j*  0,  however,  we  must  impose  on  the 
cosine  expansions  the  constraints 


u  (0)  =  v  4*  w  =0 
m  m  m 


u  ( 7r)  =  v  -  w  -  0 
m  m  m 


where 


m  ^  f ,  m 
f  even 
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w  =  v 

m  i ,  m 

i  =  odd 

Equations  (8a)  and  (8b)  will  hold  only  if  we  make  both 


v  -  0 
m 


wm  -  0  .  (10b) 

These  are  then  the  "pole  conditions"  which  the  spectral  equations  must  satisfy  for 
the*  absence  of  singularities  at  the  poles  in  Eq.  (f>).  Note  that  these  conditions  are 
consistent  with  Eqs.  (4a)  and  (4b),  the  pole  conditions  for  the  east -west  Fourier 
expansions  of  u(0,  A.). 

Just  as  Eqs.  (3a)  and  (3b)  provide  the  linkage  between  a  two-dimensional  PDE 
and  a  set  of  ODE,  Eqs.  (7a)  and  (7b)  enable  us  to  reduce  Eq.  (f>)  to  a  set  of  alge¬ 
braic  systems  in  the  Fourier  space  !m|  <  M,  0  <  (  <  L: 

(f  -  2)(f  -  l)vlm2tm  '  <2f2 3  +  4m2)u(  m  +  «  +  l)(f  +  2)uf+2j  m  * 

7  -2,  m  +  ^7,m  7+2,  m  * 

Here  ^  are  complex  Fourier  coefficients  of  f^(0).  In  the  derivation  of  Eq.  (11), 
we  have 

1.  Made  use  of  the  identity  2  sin"  0  1-  cos  20  to  row  rite  Eq.  ('>)  in  the  form 

uM  (0)  -  cos  20u"  (0)  +  sin  20u'  (0)  -  2m*'  n  (0)  (1  -  cos  2 0)f  (0)  .  (12) 

m  m  m  m  m  v  ' 


2.  Expanded  ^(0),  u^(0),  um(0)  and  f  (0)  in  Eq.  (12)  as  truncated  series 
in  the  form  of  EJqs.  (7a)  and  (7b). 

3.  Made  use  of  the  trigonometric  identities 

2  cos  20  cos  f0  -  cos  (f  -  2)0  +  cos  (t  +  2)0 
2  sin  20  sin  t  0  -  cos  (f  -  2)0  -  cos  (f  +  2)0 
2  cos  20  sin  I  0  -  sin  (f  -  2)0  +  sin  (t  +  2)0 
2  sin  20  cos  tO  -sin  (t  -  2)0  +  sin  ( t  f~  2)0  , 
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4.  Invoked  the  Galerkin  approximation  in  the  manner  that  Eq.  (5)  was  derived 
from  Eqs.  (2a)  and  (2b)  and  dropped  the  (L  +  l)th  and  (L  +  2)th  components  which 
are  outside  the  range  of  our  series  expansions  defined  by  Eqs.  (7a)  and  (7b). 

We  now  discuss  the  pole  conditions  for  Eq.  (6)  in  the  context  of  spectral  alge¬ 
braic  systems  Eq.  (11).  As  mentioned  previously,  Eqs.  (2a)  and  (2b)  prescribe  u 
on  a  sphere  only  to  within  an  additive  constant.  To  determine  u  uniquely,  we  need 
to  have  one  additional  piece  of  information  on  u.  Without  loss  of  generality,  we 
shall  impose  the  condition 

u(0  =  0)  +  u(A  -  ti)  -  constant  -  0  .  (13) 

In  Fourier  space  this  condition  becomes,  through  Eqs.  (4a)  and  (7a). 


vo  =  E  U1, 0  =  0  •  (14) 

i -even 

This  is  consistent  with  the  pole  conditions  v  for  the  cases  of  even  m  *  0  given  in 
Eqs.  (10a)  and  (10b).  In  Solving  Eq.  (11),  we  shall  use  Eq.  (14)  to  compute  q, 
the  constant  term  in  our  double  Fourier  series  expansions.  For  even  m  ^  0, 

Uq  are  computed  via  the  pole  conditions  given  in  Eq.  (10): 


u0,  m  £  uf ,  m 

£ -even 
f  *0 


For  odd  m,  un  ,  uT  are  easily  found  from  Eq.  (7b)  to  be 
0,  m*  L,  m  J 


a0,  m 


JL,  m 


(15) 


(16) 


4.  TRIDIAGONAL  SYSTEMS 

With  the  problem  of  pole  conditions  settled,  we  are  now  ready  to  solve  linear 
algebraic  systems  given  in  Eq.  (11).  We  note  first  that  for  a  given  m,  the  even  l 
and  odd  f  components  are  uncoupled.  We  may  therefore  break  the  algebraic'  sys¬ 
tem  for  each  m  into  two  independent  subsystems.  Furthermore,  written  in  matrix 
form,  these  subsystems  have  a  very  desirable  feature  in  that  their  coefficient 
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matrices  are  tridiagonal.  Very  efficient  algorithms  are  available  for  the  solution 
of  such  tridiagonal  systems  (for  example,  Varga*®).  We  write  Kq.  (11)  in  matrix 
notations: 


r  i 

CM 

3 

h2 

u4 

_ 

h4 

UL-2 

hL-2 

UL 

h, 

m  i 

r 

3 

( _ 

h:» 

uL-:i 

• 

U ! .  -  1 

hI.-l 

_  — 

1 

—  — 

-•  m 


m 


m 


(17a) 


(17b) 


where? 


!  m!  <  M  , 

n(  (t  -  2)(f  -  1) 

b,  -  2f2  -  4ni2 

f ,  m 

c’l  (I  4  l)(f  +  2) 


( ,  m 

* 

2,  m 


-  f 


t  -2. rn 


2f  -  f 
f ,  m  1  +  2.  m 


h.,  -  (-nm  rn 

2,  m  0,  m 


In  Eqs.  (17a)  and  (17b),  we  have*  dropped  the  (L  1  l,  m)  and  (1.  f  2,  rn)  compo¬ 
nents  because  they  are  beyond  the  range  of  our  series  expansions.  On  the  other 

hand,  the  components  (-2,  m)  and  (-l,m)  have  been  retained  because  u  (0)  f  (0) 

'  m  ’  ni 

are  27r-periodic  in  0.  Had  we  used  expansions  over  0  <  0  <  2*  instead  of  using  half¬ 
range  expansions  as  in  Kqs.  (7a)  and  (7b),  then  any  component.  1  f  ]  *<  1,  would  have 


10.  Varga,  R.S.  (1902)  Matrix  Iterative  Analvsis,  Prentice-Hall  Englewood  C  liffs 
N.J.,  p.  195.  ‘ 


been  *A  ithin  the  spectral  range.  We  retained  the  (-2,  m),  (-1,  m)  components  by 
tlie  use  of  the  identity 


f 


-f ,  m 


(-l)m  L 


m 


For  even  m,  we  have  E  unknowns  for  each  m,  For  odd  m,  Uj  ^  are  known  before¬ 
hand  to  be  zero,  In  this  ease,  we  simply  s€‘t  u,  0  and  discard  the  last  mem* 

l.,mf 

bet*  in  Eq.  (17a).  Note  that  the  components  Uq  m,  ;  m|  <  M,  are  not  included  in 
Eqs.  (17a)  and  (17b).  These  are  to  be  obtained  through  the  pole  conditions  dis¬ 
cussed  in  Section  3  :  For  m  uq  o  given  by  Eq.  (14), 


u0,  0  E  uf.O  : 

( 'even 
(  5*0 

for  even  m  *  0.  urt  are  given  bv  Eq.  (15).  For  odd  m,  u~  are  of  course 
*  0,  m  s  '  0,  m 

known  from  Eq.  (IB)  to  be  zero. 

So  far  we  have  demonstrated  that  the  Poisson  equation  on  a  sphere  can  be 
approximated  by  a  number  of  t  ridiagonal  systems  in  the  Fourier  space.  At.  this 
point,  it  is  appropriate  to  examine  tin*  properties  of  the  tridiagonal  systems.  A 
little  simple  arithmetic  with  ^  and  c^  reveals  that  the  coefficient  matrices 

of  these  systems  are  strictly  diagonally  dominant  for  the  cases  'ml  >  2.  Unique 
solutions  therefore  exist  for  these  cases.  For  m  0,  ±1,  the  situation  is  less 
clear  cut.  But  it  can  be  shown  that  the  determinants  of  the  coefficient  matrices 
for  these  cases  are  all  nonzero;  these  algebraic  systems  therefore  also  have 
unique  solutions. 

5.  TEST  COMPUTATIONS 


As  mentioned  earlier,  very  efficient  algorithms  are  available  for  solution  of 
tridiagonal  systems  such  as  Eqs.  (17a)  and  (17b).  We  shall  report  in  this  section 
some  results  of  our  numerical  test  computations.  Before  doing  so,  let  us  pause 
for  a  moment  to  recapitulate  our  numerical  procedure.  We  shall  also  indicate 
[in  brackelsl  the  approximate  minimum  number  of  arithmetic  operations  required 
for  each  step: 

1,  For*  each  i  1  through  J  -  1,  tr  ansform  f.  ,  to  obtain  complex  Fourier 

1 ,  K 

coefficients  f^ffl  );  set  f0(P)  f(P),  P  0  or  IT.  [ K(J  -  1)  log2  Kl 

2.  For  a  given  even  (odd)  m,  perform  cosine  (sine)  transform  on  f^(B.)  to 
obtain  f^  .  (J  log.;  J  1 


3.  Compute  in  Eqs.  (17a)  and  (17b).  [ 2 J 1 

4.  Solve  Eqs.  (17a)  and  (17b)  for  u.  ,  using  an  algorithm  for  tridiagonal 

10  f,  m 

systems.  [4  J  ] 

5.  For  a  given  even  (odd)  m,  perform  inverse  cosine  (sine)  transform  tin 

uf  to  obtain  u  (0.).  [  J  log9  J] 

J  “  Ilf 

6.  Repeat  steps  (2)  through  (5)  for  all  |m|  ^  M.  [Multiply  each  count  in 
steps  (2)  through  (5)  by  K] 

7.  Inverse  transform  j  s  lf  J  -  1,  to  obtain  u.  set  u(P)  =  Uq(P). 

[K(J  -  1)  log2  K] 

Thus,  for  a  latitude-longitude  grid  with  2J  X  (J  -  1)  +  2  data  points  (K  =  2J)f 

2 

the  approximate  number  of  arithmetic  operations  is  12J  (1  +  log£  or  a^out 

fi  (1  +  log2  J)  operations  per  data  point.  Here  we  assume  that  a^f  b^  c(  have 
been  precomputed  and  J  =  2^,  where  p  is  a  positive  integer.  An  arithmetic  opera¬ 
tion  is  defined  as  a  multiplication  followed  by  an  addition  in  the  real  domain.  These 
counts  are  comparable  to  operation  counts  for  the  Fourier  solution  of  the  Poisson 
equation  on  a  rectangle  (for  example,  Swarztrauber  *), 

The  test  function  used  for  program  checkout  is 


f(0,  X)  r  -  (m  +  l)(m  +  2)  cos  0  sinm  0  c os  m(X  -  d  )  -  2  cos  0  , 

m  =  m  ^ 

where  are  random  numbers  within  the  range  (0,  2^).  For  this  forcing  function, 
Eqs.  (2a)  and  (2b)  has  the  exact  solution 

m2 

UA(0,  *)  =  5"]  cos  0  sin111  0  cos  m(X  -  d  )  +  cos  0  . 

A  *  m 

m  •-  in  ^ 

With  the  exact  solutions  known,  a  normalized  error  norm  defined  by 

H  k  M  -  1 1 U  -  uA  I !  / 1 1  uA  1 1 

may  then  be  considered  as  a  measure  of  the  accuracy  of  our  numerical  solution  u. 
The  number  of  digits  of  accuracy  in  u  is  then  given  by 

11.  Swarztrauber,  P.  N.  (1977)  The  methods  of  cyclic  reduction,  Fourier  analysis 
and  the  FACR  algorithm  for  the  discrete  solution  of  Poisson’s  equation  on 
a  rectangle,  SIAM  Rev.  19:490-501. 


Z  =  ‘log10^  I  El  i  • 


Several  sets  of  computations  have  been  made  using  a  CDC  6600  computer  for 
various  grid  resolutions  and  various  values  of  m^  and  m2.  Some  sample  results 
are  tabulated  in  Table  1  for  cases  in  which  the  data  contain  only  long  waves 
(m^  -  1,  mg  ■  2).  We  see  that  1 1  -,  12-digit  accuracies  are  retained  in  our  numer¬ 
ical  solutions  for  the  10°  X  10°,  5°  X  5°,  2.8^  X  2.  8°  resolutions  tested.  Note  that 
excessive  resolutions  have  been  used  in  these  computations  in  the  sense  that  the 
tabulated  resolutions  are  capable  of  resolving  waves  having  wave  numbers  up  to 
18,  38,  and  84,  respectively,  while  the  data  contain  only  wave  components  up  to 
wave  number  2.  For  the  smooth  function  tested,  the  slow  degradation  of  accuracy 
with  increasing  resolution  is  due  to  the  increase  in  the  value  of  the  condition  num¬ 
ber  as  the  order  of  the  coefficient  matrices  in  Eqs.  (17a)  and  (17b)  increases.  The 
central  processor  time  required  for  each  solution  is  given  in  seconds  in  the  last 
column  in  Table  l.  It  should  be  noted  here  that  the  FFT  routine  used  is  one  that 
can  accommodate  data  sets  not  necessarily  having  2^  pieces  of  data  (p  is  a  positive 
integer).  This  special  feature  of  course  requires  an  operation  count  higher  than 
Odogg  J)  in  the  transforms  for  cases  where  J  *  2^. 

Table  2  shows  the  accuracy  of  our  numerical  results  for  a  63  X  128  grid  as  a 
function  of  mg,  which  may  be  taken  here  as  the  number  of  components  contained  in 
our  test  function.  It  is  apparent  that  there  is  also  a  gradual  degradation  in  accu¬ 
racy  as  the  function  becomes  less  smooth.  However,  even  for  a  function  containing 
components  up  to  the  Nyquist  frequency,  the  results  still  possess  more  than  9-digit 
accuracy.  Also  given  in  Table  2  is  the  normalized  error  for  an  experiment  with 
mg  =  65.  In  this  case,  one  of  the  components  has  a  wavelength  of  less  than  2 -grid 
length.  The  less  than  3 -digit  accuracy  in  the  results  is  actually  not  as  bad  as  it 
may  seem.  What  is  demonstrated  here  is  simply  that  for  any  forcing  function 
having  components  with  wavelengths  less  than  2-grid  interval,  these  components 
are  simply  truncated  by  the  transform  process.  For  components  w  ith  wavelength 
longer  than  2-grid  interval,  the  deterioration  of  accuracy  of  the  solution  with  de¬ 
creasing  wavelength  is  much  more  gradual. 
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Table  1.  Accuracy  as  a  Function  of  N  for  (in  ^ 


J 

N 

Z 

(digits) 

T 

(seconds) 

18 

9 

11.  07 

0.  59 

3  8 

18 

10.  94 

2.  35 

84 

32 

10.  88 

3.  98 

N:  the  size  of  the  tridiagonal  system 


Table  2.  Accuracy  as  a  Func  tion 
of  n^2  for  a  (83  X  128)-grid 


m2 

Z 

(digits) 

8 

10.  81 

32 

10.  52 

84 

oc 

o 

85 

2.  92  1 

_ i 

n^:  the  highest  wavenumber 
component  contained  in 
the  forcing  function 
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